Cel projektu

Problemem, który zostanie poddany analizie w prezentowanym raporcie jest wpływ różnych czynników na rozmiar śledzia oceanicznego wyławianego w Europie. Na przestrzeni ostatnich lat zauważono, że rozmiar ryby stopniowo spada. Zbiór analizowanych danych przedstawia pomiary śledzi oraz warunków, w jakich żyły na przestrzeni ostatnich 60 lat. W ramach połowu jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.

Podsumowanie

Podczas przeprowadzania analizy w pierwszym kroku uzupełnione zostały niepełne dane. Następnie w postaci histogramów zaprezentowane zostały rozkłady wartości poszczególnych atrybutów. Przedstawienie atrybutów w postaci macierzy korelacji pozwoliło zobserwować silną zależność między długością śledzi, a temperaturą przy powierzchni wody.

Podczas tworzenia modelu predykcyjnego, który na podstawie zawartych w zbiorze danych cech miał przewidziec rozmiar śledzi wykorzystany algorytm RandomForest. Model został nauczony na treningowym zbiorze danych stanowiącym 80% wszystkich posiadanych danych, a następnie przetestowany na pozostałych 20% danych.

Po dokonaniu analizy możemy zaobserwować, że najbardziej wyraźny wpływ na wielkość wyławianych śledzi ma temperatura przy powierzchni wody (sst) - jej stopniowy wzrost bezpośrednio przyczynił się do spadku średniej wielkości śledzi. Dodatkowo bardzo ważnym czynnikiem okazała się dostępność planktonu - szczególnie widłonogów, której spadek rownież wpłynął na rozmiar śledzi.

Ładowanie bibliotek

library(ggplot2)
library(DT)
library(kableExtra)
library(dplyr)
library(plotly)
library(caret)
library(gridExtra)    
library(randomForest)

Powtarzalność

Kod zapewniający powtarzalność eksperymentów

set.seed(23)

Ładowanie danych

Kod pozwalający na załadowanie danych z pliku csv. Dane zawierają 52583 rekordy opisane 16 atrybutami:

Kolumna Opis
length długość złowionego śledzia [cm]
cfin1 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]
cfin2 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]
chel1 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]
chel2 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]
lcop1 dostępność planktonu [zagęszczenie widłonogów gat. 1]
lcop2 dostępność planktonu [zagęszczenie widłonogów gat. 2]
fbar natężenie połowów w regionie [ułamek pozostawionego narybku]
recr roczny narybek [liczba śledzi]
cumf łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]
totaln łączna liczba ryb złowionych w ramach połowu [liczba śledzi]
sst temperatura przy powierzchni wody [°C]
sal poziom zasolenia wody [Knudsen ppt]
xmonth miesiąc połowu [numer miesiąca]
nao oscylacja północnoatlantycka [mb] 
raw_data <- read.csv(url("http://www.cs.put.poznan.pl/dbrzezinski/teaching/sphd/sledzie.csv"), na.strings = "?")
kable(head(raw_data, 10))%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
0 23.0 0.02778 0.27785 2.46875 NA 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
1 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
2 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
3 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
4 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
5 22.0 0.02778 0.27785 2.46875 21.43548 2.54787 NA 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
6 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
7 23.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
8 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
9 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8

Statystyki

Sekcja prezentująca podstawowe statystyki dla zbioru danych:

Atrybut Wartość
Liczba kolumn 16
Liczba wierszy 52582 
summary(raw_data[,2:16])
##      length         cfin1             cfin2             chel1       
##  Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469  
##  Median :25.5   Median : 0.1111   Median : 0.7012   Median : 5.750  
##  Mean   :25.3   Mean   : 0.4458   Mean   : 2.0248   Mean   :10.006  
##  3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936   3rd Qu.:11.500  
##  Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000  
##                 NA's   :1581      NA's   :1536      NA's   :1555    
##      chel2            lcop1              lcop2             fbar       
##  Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
##  1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
##  Median :21.673   Median :  7.0000   Median :24.859   Median :0.3320  
##  Mean   :21.221   Mean   : 12.8108   Mean   :28.419   Mean   :0.3304  
##  3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
##  Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
##  NA's   :1556     NA's   :1653       NA's   :1591                     
##       recr              cumf             totaln             sst       
##  Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77  
##  1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.60  
##  Median : 421391   Median :0.23191   Median : 539558   Median :13.86  
##  Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87  
##  3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16  
##  Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73  
##                                                        NA's   :1584   
##       sal            xmonth            nao          
##  Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Median :35.51   Median : 8.000   Median : 0.20000  
##  Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##  3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##  Max.   :35.61   Max.   :12.000   Max.   : 5.08000  
## 
no_na_data <- sum(complete.cases(raw_data))
na_data <- nrow((raw_data)) - no_na_data

Przetwarzanie brakujących danych

Analiza brakujących danych

Możemy zauważyć, że w badanym zbiorze danych aż 10094 rekordów posiada brakujące dane. Niestety usunięcie takiej ilości rekordów z badanego zbioru mogłoby wpłynąć negatywnie na jakość analizy dlatego konieczne jest uzupełnienie brakujących wartości. Analizując przedstawione statystyki możemy zaobserwować, że jedynie niektóre kolumny posiadają brakujące wartości. Są to dane dotyczące dostępności planktonu oraz temperatury przy powierzchni wody. To właśnie dane w tych cechach należy uzupełnić.

Przyglądając się danym możemy zauważyć, że braki występują w sąsiedztwie danych kompletnych. Wynika to prawdopodobnie z faktu, że sąsiadujące rekordy pochodzą z tego samego połowu, a braki w danych wynikają z nieuwagi podczas ich zbierania.

Uzupełnienie danych

Skrypt uzupełniający dane pobierany jest z pliku “ZED_grouping.R”.

Dla każdego rekordu z brakującą wartością cechy pobierane jest najbliższe 6 kompletnych rekordów. Na podstawie pobranej grupy wyliczana jest mediana brakującej wartości, która następnie wstawiana jest jako wartość badanej cechy w niekompletnym rekordzie. Czynność powtarzana jest dla każdej cechy zawierającej brakujące dane.

# Załadowanie skryptu uzupełniającego dane
if(!exists("fillMissingValues", mode="function")) source("ZED_grouping.R")

filled_data <- fillMissingValues(raw_data)
filled_data

Analiza wartości atrybutów

ggplot(filled_data, mapping = aes(x = cfin1)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = cfin2)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = chel1)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = chel2)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = lcop1)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = lcop2)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = fbar)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = recr)) + geom_histogram(binwidth = 60)

ggplot(filled_data, mapping = aes(x = cumf)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = totaln)) + geom_histogram(binwidth = 60)

ggplot(filled_data, mapping = aes(x = sst)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = sal)) + geom_histogram(binwidth = 0.005)

ggplot(filled_data, mapping = aes(x = xmonth)) + geom_histogram(binwidth = 0.01)

ggplot(filled_data, mapping = aes(x = nao)) + geom_histogram(binwidth = 0.05)

Korelacja między atrybutami

Korelacja między atrybutami została przedstawiona w postaci macierzy korelacji.

kable(cor(filled_data), digits = 2)%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))%>%
  row_spec(0,angle = -45)
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
X 1.00 -0.34 -0.15 0.06 -0.16 0.05 -0.23 0.04 0.09 0.00 0.23 -0.36 0.36 -0.06 0.00 0.41
length -0.34 1.00 0.08 0.10 0.22 -0.01 0.24 0.05 0.25 -0.01 0.01 0.10 -0.45 0.03 0.01 -0.26
cfin1 -0.15 0.08 1.00 0.15 0.09 0.20 0.12 0.21 -0.06 0.12 -0.05 0.13 0.01 0.13 0.01 0.01
cfin2 0.06 0.10 0.15 1.00 0.00 0.31 -0.04 0.65 0.15 -0.10 0.34 -0.22 -0.24 -0.08 0.01 -0.01
chel1 -0.16 0.22 0.09 0.00 1.00 0.29 0.95 0.25 0.16 -0.05 0.07 0.17 -0.22 -0.15 0.05 -0.51
chel2 0.05 -0.01 0.20 0.31 0.29 1.00 0.17 0.88 0.03 0.00 0.26 -0.38 0.01 -0.22 0.07 -0.06
lcop1 -0.23 0.24 0.12 -0.04 0.95 0.17 1.00 0.15 0.10 0.00 -0.01 0.27 -0.26 -0.10 0.03 -0.55
lcop2 0.04 0.05 0.21 0.65 0.25 0.88 0.15 1.00 0.05 0.00 0.29 -0.30 -0.12 -0.18 0.06 -0.04
fbar 0.09 0.25 -0.06 0.15 0.16 0.03 0.10 0.05 1.00 -0.24 0.82 -0.51 -0.18 0.04 0.01 0.07
recr 0.00 -0.01 0.12 -0.10 -0.05 0.00 0.00 0.00 -0.24 1.00 -0.26 0.37 -0.20 0.28 0.02 0.09
cumf 0.23 0.01 -0.05 0.34 0.07 0.26 -0.01 0.29 0.82 -0.26 1.00 -0.71 0.03 -0.10 0.03 0.23
totaln -0.36 0.10 0.13 -0.22 0.17 -0.38 0.27 -0.30 -0.51 0.37 -0.71 1.00 -0.29 0.15 -0.03 -0.39
sst 0.36 -0.45 0.01 -0.24 -0.22 0.01 -0.26 -0.12 -0.18 -0.20 0.03 -0.29 1.00 0.01 -0.01 0.51
sal -0.06 0.03 0.13 -0.08 -0.15 -0.22 -0.10 -0.18 0.04 0.28 -0.10 0.15 0.01 1.00 -0.02 0.13
xmonth 0.00 0.01 0.01 0.01 0.05 0.07 0.03 0.06 0.01 0.02 0.03 -0.03 -0.01 -0.02 1.00 0.00
nao 0.41 -0.26 0.01 -0.01 -0.51 -0.06 -0.55 -0.04 0.07 0.09 0.23 -0.39 0.51 0.13 0.00 1.00

Analiza rozmiaru śledzi w czasie

Kod umożliwiający analizę rozmiaru śledzi w czasie.

W celu analizy rozmiaru rozmiaru śledzi w czasie, na podstawie informacji o rocznym narybku oraz wiedzy, że dane występują w porządku chronologicznym utworzona została dodatkowa kolumna “year”.

if(!exists("add_year", mode="function")) source("ZED_grouping.R")

mean_sample <- filled_data[,c("recr","length")]
recr_groups <- aggregate(mean_sample, by = list(mean_sample$recr), FUN = mean)
year_data <- add_year(filled_data, recr_groups)
aggreagated_years <- aggregate(year_data, by= list(year_data$year), FUN = mean)%>%select(year, length, lcop2, lcop1, chel2, sst)

interactive_plot <- ggplot(aggreagated_years, mapping = aes(x = year, y = length))+geom_smooth()
ggplotly(interactive_plot)

Regresor

W następnym kroku podjęta została próba utworzenie regresora, który przewidzi rozmiar śledzi. Dane zostały podzielone na zbiór uczący oraz testowy w proporcji 80%/20%. Wykorzystany został algorytm RandomForest optymalizujący model według miary RMSE.

sig_attr_selection <- filled_data%>%select(length, cfin1, cfin2, chel1, chel2, lcop1, lcop2, fbar, recr, cumf, totaln, sst, sal, nao)
inTraining <- createDataPartition(y = sig_attr_selection$length, p=.8, list=FALSE)

training_data <- sig_attr_selection[inTraining,]
test_data <- sig_attr_selection[-inTraining,]

rfGrid <- expand.grid(mtry = 1:10)
ctrl <- trainControl(
    method = "repeatedcv",
    number = 5,
    repeats = 5)

fit <- train(length ~ .,
             data = training_data,
             method = 'rf',
             trControl = ctrl,
             metric = "RMSE",
             tuneGrid=rfGrid,
             importance = TRUE,
             ntree=10)

ggplot(fit) + theme_bw()

rfClasses <-predict(fit, newdata = test_data[-1])
rmse <- RMSE(pred = rfClasses, obs = test_data$length)
r2 <- R2(pred = rfClasses, obs = test_data$length)

Algorytm osiągnął najlepszy rezultat dla parametru “mtry” równego 3.

Miara Wartość
RMSE 1.1825109
R2 0.4916135

Analiza ważności atrybutów w wyuczonym modelu

Na podstawie wyuczonego modelu oszacowane zostały wagi każdego z atrybutów zbioru danych.

attr_importance <- varImp(fit, scale = FALSE)
ggplot(attr_importance)

Możemy przyjrzeć się zmiane najbardziej istotnych cech w czasie.

lcop2_plot <- ggplot(filled_data, mapping = aes(x = X, y = lcop2))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie widłonogów gat. 2]")
ggplotly(lcop2_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
chel2_plot <- ggplot(filled_data, mapping = aes(x = X, y = chel2))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]")
ggplotly(chel2_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
lcop1_plot <- ggplot(filled_data, mapping = aes(x = X, y = lcop1))+geom_smooth()+ggtitle("Dostępność planktonu [zagęszczenie widłonogów gat. 1]")
ggplotly(lcop1_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
fbar_plot <- ggplot(filled_data, mapping = aes(x = X, y = fbar))+geom_smooth()+ggtitle("Ułamek pozostawionego narybku")
ggplotly(fbar_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
sst_plot <- ggplot(filled_data, mapping = aes(x = X, y = sst))+geom_smooth()+ggtitle("Temperatura przy powierzchni wody [°C]")
ggplotly(sst_plot)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Oczywiście wpływ na wielkość wyławianych śledzi miało wiele czynników, jednak wśród najbardziej wpływowych wyróżnić możemy dostępność planktonu oraz temperaturę przy powierzchni wody.